************************************************************************************
* 'Gender Dysphoria and Psychological Functioning in Adolescents Treated with GnRHa: 
*  Comparing Dutch and English Prospective Studies'
* Archives of Sexual Behavior, accepted 5 June 2020
* Michael Biggs, Department of Sociology and St Cross College, University of Oxford
************************************************************************************

version 16

local msize		"medlarge"
local fontsize  "12pt"
local Englishgrey "gs9"

cd "/Users/michael/Documents/Gender critical/NHS/Puberty blockers/analysis"
import delimited using rawdata.csv, clear varnames(1) encoding("utf-8") case(preserve)
foreach v of varlist t0* t1* {
	replace `v' = round(`v', .01)
	}
format t0* t1* %4.2f

gen order = .
replace order = 11 if measure=="CBCL-T"
replace order = 10 if measure=="CBCL-I"
replace order = 9  if measure=="CBCL-E"
replace order = 8  if measure=="YSR-T"
replace order = 7  if measure=="YSR-I"
replace order = 6  if measure=="YSR-E"
replace order = 5  if measure=="CGAS"
replace order = 4  if measure=="UGDS"
replace order = 3  if measure=="BIS-1"
replace order = 2  if measure=="BIS-2"
replace order = 1  if measure=="BIS-0"

gen n = .
replace n = F_errdf +1 if sample=="English"
replace n = round(total_n * 37/70, 1) if sample=="Dutch" & sex=="F" 
replace n = round(total_n * 43/70, 1) if sample=="Dutch" & sex=="M" 

by sex sample, sort: summ t0_mean if measure=="UGDS" 
foreach var in t0_mean t0_sd t1_mean {
	replace `var' = `var' / 12 if measure=="UGDS" & sample=="Dutch"
	}
by sex sample, sort: summ t0_mean if measure=="UGDS" 

******************************************************************************
******************************************************************************

*****************************
* COMPARE SAMPLES AT BASELINE 
*****************************

gen t0_lb = .
gen t0_ub = .
levelsof measure, local(mlevels)
foreach m of local mlevels {
	foreach s in M F {
		foreach c in Dutch English {
			foreach x in n t0_mean t0_sd  {
*				di "`x' `c' `s' `m'"
				quietly summ `x' if sample=="`c'" & sex=="`s'" & measure=="`m'"
				local `x'_`c' = r(mean)
				}
			quietly cii means `n_`c'' `t0_mean_`c'' `t0_sd_`c''
			quietly replace t0_lb = r(lb) if sample=="`c'" & sex=="`s'" & measure=="`m'"
			quietly replace t0_ub = r(ub) if sample=="`c'" & sex=="`s'" & measure=="`m'"			
			}
		di _n "`m'--`s':" 
		quietly sdtesti `n_Dutch' `t0_mean_Dutch' `t0_sd_Dutch' `n_English' `t0_mean_English' `t0_sd_English' 
		if r(p) < .05 {
			di "  (variance differs: p = " %6.5f (r(p))
			}
		local difference = 0
		quietly ttesti `n_Dutch' `t0_mean_Dutch' `t0_sd_Dutch' `n_English' `t0_mean_English' `t0_sd_English', unequal
		if r(p) < .05 {
			di "     p (unequal) = " %6.5f (r(p)) 
			local difference = `t0_mean_Dutch' - `t0_mean_English'
			}
		if `difference' != 0	di "     Difference: " %5.3f (`difference')

		}	
	}

***************************
* FIGURES 1 AND 2: BASELINE 
***************************

gen row =  order + 6 + ((sex=="F") *8) + ((sample=="Dutch") *.25) if type=="PF"
gsort - row sample
list measure order sex sample row  if type=="PF", clean 

forvalues r = 11/26 { 
	if `r'==18 		{
		local label`r' `""{bf:Male}""'
		}
	else if `r'==26 {
		local label`r' `""{bf:Female}""'
		}
	else {
		quietly levelsof measure_long if row==`r', local(label`r')
		}
	local PFylabels "`PFylabels' `r' `label`r'' " 
	}
di `"`PFylabels'"'

replace row =  order + ((sex=="F") *5) + ((sample=="Dutch") *.25) if type=="GD"
gsort - row sample
list measure order sex sample row  if type=="GD", clean 

forvalues r = 1/10 { 
	if `r'==5 		{
		local label`r' `""{bf:Male}""'
		}
	else if `r'==10 {
		local label`r' `""{bf:Female}""'
		}
	else {
		quietly levelsof measure_long if row==`r', local(label`r')
		}
	local GDylabels "`GDylabels' `r' `label`r'' " 
	}
di `"`GDylabels'"'
	
replace row = row - .125

twoway ///
	(rcap t0_lb t0_ub row if type=="PF" & sample=="Dutch", horizontal lcolor(black) lwidth(medthick) msize(tiny)) ///
	(rcap t0_lb t0_ub row if type=="PF" & sample=="English", horizontal color(`Englishgrey') lwidth(medthick) msize(tiny)) ///
	(scatter row t0_mean if type=="PF" & sample=="Dutch", mcolor(black) msize(`msize')) ///
	(scatter row t0_mean if type=="PF" & sample=="English", mcolor(`Englishgrey') msize(`msize')) ///
	, ///
	ytitle("") yscale(noline range(10.5 26)) ///
	ylabel( `PFylabels', angle(0) noticks labgap(2) nogrid labsize(`fontsize')) ///
	xtitle("")   ///
	xlabel(50(10)70, grid labsize(`fontsize')) ///
	title("{bf:Fig. 1}  Psychological functioning at baseline", size(14pt) span color(black) margin(b+3) position(12)) ///
	plotregion(margin(zero)) graphregion(fcolor(white) lcolor(white) margin(l-4 b-2 t-2)) /// margin(l+0 r+0 b-1 t-1
	legend(region(lcolor(white)) position(7) order(3 4) label(3 "Dutch") label(4 "English") size(`fontsize')) ///
	ysize(6) ///
	name(Fig1_PF, replace)
graph export Fig1.jpg, as(jpg) replace

twoway ///
	(rcap t0_lb t0_ub row if type=="GD" & sample=="Dutch", horizontal lcolor(black) lwidth(medthick) msize(tiny)) ///
	(rcap t0_lb t0_ub row if type=="GD" & sample=="English", horizontal color(`Englishgrey') lwidth(medthick) msize(tiny)) ///
	(scatter row t0_mean if type=="GD" & sample=="Dutch", mcolor(black) msize(`msize')) ///
	(scatter row t0_mean if type=="GD" & sample=="English", mcolor(`Englishgrey') msize(`msize')) ///
	, ///
	ytitle("") yscale(noline range(0.5 10)) ///
	ylabel(`GDylabels', angle(0) noticks labgap(2) nogrid labsize(`fontsize')) ///
	xtitle("")   ///
	xlabel(, labsize(`fontsize') grid) ///
	title("{bf:Fig. 2}  Gender dysphoria at baseline", size(14pt) span color(black) margin(b+3) position(12)) ///
	plotregion(margin(zero)) graphregion(fcolor(white) lcolor(white) margin(l+5 b-2 t-2)) /// margin(l+0 r+0 b-1 t-1
	legend(region(lcolor(white)) position(7) order(3 4) label(3 "Dutch") label(4 "English") size(`fontsize')) ///
	ysize(4)  ///
	name(Fig2_GD, replace)
graph export Fig2.jpg, as(jpg) replace

******************************************************************************
******************************************************************************

*********************
* COMPUTE IMPROVEMENT 
*********************

gen rawchange = (t1_mean - t0_mean) // / t0_sd
gen change = rawchange / t0_sd
format change rawchange %4.2f

gen rawchange_sd = abs(rawchange) * sqrt(n) / sqrt(F) if sample=="English"
gen F_new = invFtail(1,n-1,p) if F<0
list measure sex F F_new if F<0
replace rawchange_sd = abs(rawchange) * sqrt(n) / sqrt(F_new) if F<0

gen new_p = .
gen change_lb = .
gen change_ub = .
local max = _N
forvalues i = 1/`max' {
	if rawchange_sd[`i'] !=. {
		local n = n[`i']
		local change = rawchange[`i']
		local change_sd = rawchange_sd[`i']
		quietly ttesti `n' `change' `change_sd' 0   
		quietly replace new_p = r(p) if _n==`i'
		quietly cii means `n' `change' `change_sd'
		quietly replace change_lb = r(lb)/t0_sd if _n==`i'
		quietly replace change_ub = r(ub)/t0_sd if _n==`i'			
		}
	}
twoway scatter p new_p if new_p!=., mlabel(measure)
list measure t0_sd rawchange rawchange_sd if F<0

twoway scatter rawchange_sd t0_sd if new_p!=., mlabel(measure)

di normal(.5)
di normal(1)
di normal(-.75)

* Transform so that positive is now improvement
foreach v of varlist change* { 
	replace `v' = - `v' if measure!="CGAS"  
	}

list measure sex new_p if new_p<.06, clean noobs


by measure sample, sort: egen change_F = total(cond(sex=="F", change, .))
by measure sample, sort: egen change_M = total(cond(sex=="M", change, .))
replace change_F = . if sex=="M"
replace change_M = . if sex=="M"
list sex change* if measure=="UGDS" & sample=="Dutch"
correl change_F change_M if sample=="Dutch"
correl change_F change_M if sample=="English"

by measure sex, sort: egen change_D = total(cond(sample=="Dutch", change, .))
by measure sex, sort: egen change_E = total(cond(sample=="English", change, .))
replace change_D = . if sample!="English"
replace change_E = . if sample!="English"
list sample change* if measure=="UGDS" & sex=="M"
correl change_D change_E if sex=="F"
correl change_D change_E if sex=="M"
correl change_D change_E if type=="GD"
correl change_D change_E if type=="PF"

******************************
* FIGURES 3 AND 4: IMPROVEMENT 
******************************

replace row = row + .1 if sample=="Dutch"
replace row = row - .1 if sample=="English"
twoway ///
	(bar change row if type=="PF" & sample=="Dutch", horizontal color(black) barwidth(.45) ) ///
	(bar change row if type=="PF" & sample=="English", horizontal color(`Englishgrey') barwidth(.45)   ) /// 
	(rcap change_lb change_ub row if type=="PF" & sample=="English", horizontal lcolor(gs7) msize(tiny) lwidth(medthick)) /// 
	, ///
	ytitle("") yscale(noline range(10.5 26)) ///
	ylabel(`PFylabels', angle(0) noticks labgap(2) nogrid labsize(`fontsize')) ///
	xtitle("Standard deviation improvement", size(`fontsize'))   ///
	xlabel(-1 "-1" 0 "0" 1 "1", labsize(`fontsize') grid) ///
	title("{bf:Fig. 3}  Change in psychological functioning after treatment", size(14pt) span color(black) margin(b+3) ) ///
 	plotregion(margin(zero)) graphregion(fcolor(white) lcolor(white) margin(l+10 b-2 t-2)) /// margin(l+0 r+0 b-1 t-1
	legend(region(lcolor(white)) position(7) order(1 2) label(1 "Dutch") label(2 "English") size(`fontsize')) ///
	ysize(6) ///
	name(Fig3_PF, replace)
graph export Fig3.jpg, as(jpg) replace

	
twoway ///
	(bar change row if type=="GD" & sample=="Dutch", horizontal color(black) barwidth(.45) ) ///
	(bar change row if type=="GD" & sample=="English", horizontal color(`Englishgrey') barwidth(.45)   ) /// lwidth(vvthick)
	(rcap change_lb change_ub row if type=="GD" & sample=="English", horizontal lcolor(gs7) msize(tiny) lwidth(medthick)) ///
	, ///
	ytitle("") yscale(noline range(0.5 10)) ///
	ylabel(`GDylabels', angle(0) noticks labgap(2) nogrid labsize(`fontsize')) ///
	xtitle("Standard deviation improvement", size(`fontsize'))   ///
	xlabel(-2 "-2" -1 "-1" 0 "0" 1 "1", labsize(`fontsize') grid) ///
	title("{bf:Fig. 4}  Change in gender dysphoria after treatment", size(14pt) span color(black) margin(b+3)) ///
	plotregion(margin(zero)) graphregion(fcolor(white) lcolor(white) margin(l+10 b-2 t-2)) /// margin(l+0 r+0 b-1 t-1
	legend(region(lcolor(white)) position(7) order(1 2) label(1 "Dutch") label(2 "English") size(`fontsize')) ///
	ysize(4) ///
	name(Fig4_GD, replace)
graph export Fig4.jpg, as(jpg) replace

******************************************************************************
******************************************************************************

********************
* DESCRIPTIVE TABLES 
********************

* TABLE 1

gsort - type + order sample sex
list sample sex n total_n in 1/8, clean
capture: file close t
file open t using "../Latex/table1.tex", write text replace
file write t "  &        Dutch  & English  & English  & English      \\" _n
file write t "Measure  & total  & total    & females  & males \\" _n
file write t "\midrule" _n

count
local N = r(N)
forvalues i = 1(4)`N' {
	file write t (  subinstr( subinstr(measure_long[`i'],"—","---", .) ), "’", "'", .  ) "&"
	file write t (total_n[`i']) "&" 
	file write t (n[`i'+2]+n[`i'+3])  "&"  (n[`i'+2]) "&" (n[`i'+3])   "\\" _n
	}
file close t


* TABLE S1

file open t using "../Latex/table2.tex", write text replace
gsort - order + sample sex
count
local N = r(N)
forvalues i = 1/`N' {
	file write t (measure[`i']) "&"
	file write t (sample[`i']) "&"
	file write t (sex[`i']) "&" 
	file write t (n[`i'])  "&" // (cond(sample[`i']=="Dutch", "*", "")) 
	file write t %3.2f (t0_mean[`i']) "&" 
	file write t %3.2f (t0_sd[`i']) "&" 
	file write t %3.2f (rawchange[`i']) "&" 
	file write t %3.2f (rawchange_sd[`i']) "\\" _n 
	}
file close t
